International Cherry Blossom Prediction Competition

Exploratory Data Analysis 1

Author

Edimer David Jaramillo

Published

February 20, 2024

Document description

En este documento se lleva a cabo análisis exploratorio de datos con las variables bioclimáticas, de suelo, coberturas terrestres y huella humana. Recordemos que las variables bioclimáticas son un consolidado de los años 1970 a 2000. Este documento tiene como finalidad comparar las fechas de floración entre países, validar relaciones que han sido reportadas en literatura (por ejemplo, la relación de temperatura con DOY), determinar la existencia de efectos de interacción y establecer cuáles de las variables predictoras podrían ser eventualmente candidatas para incluir en la modelación (selección de predictores). Podríamos decir que las variables analizadas en este documento son “estáticas” o “fijas”, ya que no tenemos temporalidad de ellas, son datos consolidados en ventanas de tiempo preestablecidas.
En el siguiente documento (06-EDA2.qmd) se muestran resultados de análisis exploratorio para variables “dinámicas”, me refiero a ellas como dinámicas ya que las variables climatológicas y de fotoperíodo tienen temporalidad (datos diarios) y obtuve información desde el 01 de enero de 1981 hasta el 31 de diciembre de 2023.

Libraries and setup

Code

library(tidyverse)
library(arrow)
library(splines)
library(corrr)
library(scales)
library(FactoMineR)
library(factoextra)
library(glue)

theme_set(theme_bw() + theme(legend.position = "top"))

Data

Extraigo una variable de nombre country que tiene el país, esto lo hago con la finalidad de evaluar patrones de comportamiento por diferentes países y porque podría ser de ayuda para extrapolar a regiones donde no tenemos registros históricos de floración. Para diferenciar los datos de NPN para USA, los etiqueto en la variable country como USA-NPN y a los de Washington DC como USA-WDC.
En varios de los resultados de ese documento no se incluye USA-WDC porque sólo es una coordenada (lat=38.88535, long=-77.03863) y tampoco se muestran los resultados de USA-NPN, la razón es que las variables bioclimáticas no están para las coordenadas de USA-NPN, sin embargo, como información de clima fue suministrada por NPN, uso esta información para análisis exploratorio al final de este documento y también en el próximo (06-EDA2.qmd).

Code

df_full <- read_parquet("../external-data/df_full.parquet") |> 
  mutate(
    country = 
      case_when(
        str_detect(location, "Japan") ~ "Japan",
        str_detect(location, "kyoto") ~ "Japan",
        str_detect(location, "liestal") ~ "Switzerland",
        str_detect(location, "Switzerland") ~ "Switzerland",
        str_detect(location, "South Korea") ~ "South Korea",
        str_detect(location, "vancouver") ~ "Canada",
        str_detect(location, "washingtondc") ~ "USA-WDC",
        str_detect(location, "NPN") ~ "USA-NPN"
      )
  ) |> 
  relocate(location, country, everything())

DOY distribution by country

Code

df_full |> 
  filter(location != "NPN") |> 
  ggplot(aes(x = bloom_doy)) +
  facet_wrap(~country, ncol = 1, scales = "free_y") +
  geom_histogram(color = "black") +
  geom_rug() +
  labs(x = "Bloom DOY", y = "Count",
       title = "DOY distribution by country",
       subtitle = "Original scale")
df_full |> 
  filter(location != "NPN") |> 
  ggplot(aes(x = bloom_doy)) +
  facet_wrap(~country, ncol = 1, scales = "free_y") +
  geom_histogram(color = "black") +
  geom_rug() +
  scale_x_log10() +
  labs(x = "Bloom DOY", y = "Count",
       title = "DOY distribution by country",
       subtitle = "Logarithmic scale")

Scatter plot of all vs DOY

Los siguientes gráficos permiten evidenciar relaciones (lineales y no lineales) entre algunas de las variables bioclimáticas y el día de la floración, también se observan asociaciones entre composición de suelo, por ejemplo, pH y DOY, no son dicientes o claras las relaciones entre las variables de cobertura terreste, huella humana, pastizales, arbutos, agua y la fecha de la floración.
Con estos gráficos se valida o confirma que la temperatura (bio1, bio5, bio8, bio9, bio10 y bio11) y el día de la floración exhiben una relación que podríamos sugerir está entre lineal y cuadrática, independiente del país, a mayor temperatura se han observado floraciones más precoces, este comportamiento es similar en Japón y Suiza, no obstante, Corea del Sur presenta algunas diferencias, por ejemplo, en la variable bio11 para Japón y Suiza, la relación de tipo cuadrático es evidente, pero en Corea del sur no lo es.
Las precipitaciones (bio12 a bio19) también muestran relaciones interesantes con la fecha de floración, con ligeras diferencias entre países, por ejemplo, la precipitación anual (bio12) muestra una relación lineal inversa con el DOY, en Suiza por el contrario observamos una relación cuadrática positiva y en Corea del Sur no es clara la relación. Este resultado puede sugerir que temperatura y precipitación son factores que interactuan y tener un efecto conjunto; más adelante se exploran interacciones entre variables para determinar covariaciones de interés.
La estacionalidad de la temperatura (bio4) muestra un patrón de comportamiento diferente entre países, en Japón la relación es positiva y lineal, en Corea del Sur no es clara la asociación y en Suiza muesra una relación inversa cuadrática. Recordemos que las variables bioclimáticas (bio_) recogen información de los años 1970 a 2000, con una resolución de $1 km^2$, este hallazgo será comparado más adelante con la información diaria de temperatura para validar las relaciones que se evidencian en este gráfico.
En cuanto a la composición del suelo, el nitrógeno, el pH del agua del suelo, el contenido de carbono orgánico del suelo (SOC), la densidad aparente (BDOD) y la capacidad de intercambio catiónico (CEC) muestran tendencias en la relación con el DOY. Las otras variables parecen no tener ninguna asociación (lineal o no lineal) con la variable objetivo.
Finalmente, la información de tierras de cultivo (CROPLAND), huella humana (FOOTPRINT) y coberturas terrestres, no exhiben ningún tipo de relación que pudiera ser considerada de interés para análisis posteriores, además, algunas de estas variables (por ejemplo SHRUBS) carecen de información para las coordenadas bajo análisis.

Code

df_full |>
  filter(country == "Japan") |> 
  select(-c(location, country,  year, bloom_date)) |>
  pivot_longer(cols = -bloom_doy) |>
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point(size = 0.5, alpha = 0.25) +
  stat_bin_2d(aes(fill = ..density..),
                  geom = "raster", 
                  contour = FALSE,
                  show.legend = FALSE) +
  scale_fill_distiller(palette = 4, direction = -1) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "firebrick3",
              se = FALSE) +
  scale_x_log10() +
  labs(y = "DOY", x = "", title = "Japan")

Code

df_full |>
  filter(country == "Switzerland") |> 
  select(-c(location, country,  year, bloom_date)) |>
  pivot_longer(cols = -bloom_doy) |>
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point(size = 0.5, alpha = 0.25) +
  stat_bin_2d(aes(fill = ..density..),
                  geom = "raster", 
                  contour = FALSE,
                  show.legend = FALSE) +
  scale_fill_distiller(palette = 4, direction = -1) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "firebrick3",
              se = FALSE) +
  scale_x_log10() +
  labs(y = "DOY", x = "", title = "Switzerland")

Code

df_full |>
  filter(country == "South Korea") |>
  select(-c(location, country,  year, bloom_date)) |>
  pivot_longer(cols = -bloom_doy) |>
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point(size = 0.5, alpha = 0.25) +
  stat_bin_2d(aes(fill = ..density..),
                  geom = "raster", 
                  contour = FALSE,
                  show.legend = FALSE) +
  scale_fill_distiller(palette = 4, direction = -1) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "firebrick3",
              se = FALSE) +
  scale_x_log10() +
  labs(y = "DOY", x = "", title = "South Korea")

Code

df_full |>
  filter(location != "NPN") |>
  select(-c(location, country,  year, bloom_date)) |>
  pivot_longer(cols = -bloom_doy) |>
  ggplot(aes(x = value, y = bloom_doy)) +
  facet_wrap( ~ name, scales = "free", ncol = 4) +
  geom_point(size = 0.5, alpha = 0.25) +
  stat_bin_2d(aes(fill = ..density..),
                  geom = "raster", 
                  contour = FALSE,
                  show.legend = FALSE) +
  scale_fill_distiller(palette = 4, direction = -1) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "firebrick3",
              se = FALSE) +
  scale_x_log10() +
  labs(y = "DOY", x = "", title = "All the countries")

Non-parametric correlations

Estos perfiles de correlación muestran para cada país y todos los países, de mayor a menor correlación (Spearman). En Japón y Corea del Sur, los cerezos ubicados más al norte parecen tardarse más en la floración (colores azules), la variable latitud es la primera en el perfil de correlaciones de estos países, sin embargo, en Suiza no es la latitud el factor lineal de mayor influcencia, es la altitud seguido de variables de composición de suelo y en el top 10 están tres variables que recogen información de precipitación. Si observamos la parte baja del gráfico (colores rojos) en todos los países coincide la variable temperatura anual (bio1) en ser el factor lineal inverso de mayor magnitud, seguido de otras variables indicadoras también de temperaturas. Este gráfico es de mucha utilidad no sólo para representar relaciones de tipo lineal sino también para obtener una caracterización rápida de la relación del clima con la fenología en cerezos ubicados en diferentes países, podemos decir que si pensamos en los factores que se asociación con tardanza en la floración, los perfiles entre países son muy diferentes, no obstante, si pensamos en los factores que se asociación con la precocidad en la floración, los perfiles entre países son muy parecidos. En conclusión, florecer más rápido parece tener unas “causas” comunes y globales ligadas a la temperatura, sin embargo, el retraso en la floración podría ser más difícil de establecer y generalizar.

Code

df_full |>
  filter(country == "Japan") |>
  select(where(is.numeric), -c(shrubs)) |>
  correlate(method = "spearman") |>
  filter(term == "bloom_doy") |>
  pivot_longer(cols = -term) |>
  filter(!is.na(value)) |>
  mutate(
    name = str_replace_all(
      name,
      "wildareas-v3-1993-human-footprint",
      "human-footprint93"
    ),
    name = str_replace_all(
      name,
      "wildareas-v3-2009-human-footprint",
      "human-footprint09"
    )
  ) |>
  ggplot(aes(
    x = reorder(name, value),
    y = term,
    fill = value
  )) +
  geom_tile() +
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  theme(axis.text.x = element_text(angle = 35, hjust = 1)) +
  labs(x = "",
       y = "",
       fill = "",
       title = "Japan") +
  coord_flip()
df_full |> 
  filter(country == "Switzerland") |> 
  select(where(is.numeric), -c(shrubs)) |> 
  correlate(method = "spearman") |> 
  filter(term == "bloom_doy") |> 
  pivot_longer(cols = -term) |> 
  filter(!is.na(value)) |> 
  mutate(name = str_replace_all(name,
                                "wildareas-v3-1993-human-footprint",
                                "human-footprint93"),
         name = str_replace_all(name,
                                "wildareas-v3-2009-human-footprint",
                                "human-footprint09")) |> 
  ggplot(aes(x = reorder(name, value), y = term, fill = value)) +
  geom_tile() +
  geom_tile() +
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  theme(axis.text.x = element_text(angle = 35, hjust = 1)) +
  labs(x = "", y = "", fill = "", title = "Switzerland") +
  coord_flip()
df_full |> 
  filter(country == "South Korea") |> 
  select(where(is.numeric), -c(shrubs)) |> 
  correlate(method = "spearman") |> 
  filter(term == "bloom_doy") |> 
  pivot_longer(cols = -term) |> 
  filter(!is.na(value)) |> 
  mutate(name = str_replace_all(name,
                                "wildareas-v3-1993-human-footprint",
                                "human-footprint93"),
         name = str_replace_all(name,
                                "wildareas-v3-2009-human-footprint",
                                "human-footprint09")) |> 
  ggplot(aes(x = reorder(name, value), y = term, fill = value)) +
  geom_tile()  +
  geom_tile() +
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  theme(axis.text.x = element_text(angle = 35, hjust = 1)) +
  labs(x = "", y = "", fill = "", title = "South Korea") +
  coord_flip()
df_full |> 
  filter(location != "NPN") |> 
  select(where(is.numeric), -c(shrubs)) |> 
  correlate(method = "spearman") |> 
  filter(term == "bloom_doy") |> 
  pivot_longer(cols = -term) |> 
  filter(!is.na(value)) |> 
  mutate(name = str_replace_all(name,
                                "wildareas-v3-1993-human-footprint",
                                "human-footprint93"),
         name = str_replace_all(name,
                                "wildareas-v3-2009-human-footprint",
                                "human-footprint09")) |> 
  ggplot(aes(x = reorder(name, value), y = term, fill = value)) +
  geom_tile()  +
  geom_tile() +
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  theme(axis.text.x = element_text(angle = 35, hjust = 1),
        legend.key.size = unit(0.85, "cm")) +
  labs(x = "", y = "", fill = "", title = "All the countries") +
  coord_flip()

Interaction of variables

Con los siguientes gráficos pretendo reconocer si existen comportamientos de interacción, es decir, intento confirmar a partir de representaciones visuales si hay efectos conjuntos entre las variables predictoras. En vista de que el conjunto de variables predictoras es “grande” $(P = 41)$ generar todas las interacciones (dobles, triples o de mayor jerarquía) podría ser dispendioso y demorado. Las siguientes son las dos aproximaciones que implemento:
- 1. Inicialmente voy a graficar diagramas de dispersión con las variables que mostraron tendencias en los gráficos de dispersión y los perfiles de correlación; el color de los puntos está determinado por la variable respuesta bloom_doy.
- 2. Grafico la interacción doble para algunas variables de interés vs la variable respuesta.

Scatter plot

En estos gráficos aplico la función distinct() por latitud y longitud porque la misma coordenada (así sean en años diferentes) tiene el mismo dato de bioclimáticas, suelo y coberturas, ya que esta información no posee temporalidad. Por esa razón resumo con la mediana del bloom_doy (ya que esta variable sí cambia con el tiempo).

Code

df_full |>
  filter(country == "Japan") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = bio12, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Annual precipitation (mm)",
       color = "DOY (median)",
       title = "Japan") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(country == "Switzerland") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = bio12, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Annual precipitation (mm)",
       color = "DOY (median)",
       title = "Switzerland") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(country == "South Korea") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = bio12, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Annual precipitation (mm)",
       color = "DOY (median)",
       title = "South Korea") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)

Code

df_full |>
  filter(country == "Japan") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Nitrogen (cg/kg)",
       color = "DOY (median)",
       title = "Japan") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(country == "Switzerland") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Nitrogen (cg/kg)",
       color = "DOY (median)",
       title = "Switzerland") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))  +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(country == "South Korea") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Nitrogen (cg/kg)",
       color = "DOY (median)",
       title = "South Korea") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)

Code

df_full |>
  filter(`soc_0-5cm_mean` > 0) |> 
  filter(country == "Japan") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `soc_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Soil organic carbon (dg/kg)",
       color = "DOY (median)",
       title = "Japan") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(`soc_0-5cm_mean` > 0) |> 
  filter(country == "Switzerland") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `soc_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Soil organic carbon (dg/kg)",
       color = "DOY (median)",
       title = "Switzerland") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))  +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(`soc_0-5cm_mean` > 0) |> 
  filter(country == "South Korea") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `soc_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Soil organic carbon (dg/kg)",
       color = "DOY (median)",
       title = "South Korea") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)

Code

df_full |>
  filter(`bdod_0-5cm_mean` > 0) |> 
  filter(country == "Japan") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `bdod_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Bulk density (cg/cm³)",
       color = "DOY (median)",
       title = "Japan") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(`bdod_0-5cm_mean` > 0) |> 
  filter(country == "Switzerland") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `bdod_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Bulk density (cg/cm³)",
       color = "DOY (median)",
       title = "Switzerland") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))  +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(`bdod_0-5cm_mean` > 0) |> 
  filter(country == "South Korea") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio1, y = `bdod_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual mean temperature (°C)",
       y = "Bulk density (cg/cm³)",
       color = "DOY (median)",
       title = "South Korea") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)

Code

df_full |>
  filter(`nitrogen_0-5cm_mean` > 0) |> 
  filter(country == "Japan") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio12, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual precipitation (mm)",
       y = "Nitrogen (cg/kg)",
       color = "DOY (median)",
       title = "Japan") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(`nitrogen_0-5cm_mean` > 0) |> 
  filter(country == "Switzerland") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio12, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual precipitation (mm)",
       y = "Nitrogen (cg/kg)",
       color = "DOY (median)",
       title = "Switzerland") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))  +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(`nitrogen_0-5cm_mean` > 0) |>  
  filter(country == "South Korea") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = bio12, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(x = "Annual precipitation (mm)",
       y = "Nitrogen (cg/kg)",
       color = "DOY (median)",
       title = "South Korea") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)

Code

df_full |>
  filter(country == "Japan") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = alt, y = bio12, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(y = "Annual precipitation (mm)",
       x = "Altitude",
       color = "DOY (median)",
       title = "Japan") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(country == "Switzerland") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = alt, y = bio12, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(y = "Annual precipitation (mm)",
       x = "Altitude)",
       color = "DOY (median)",
       title = "Switzerland") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))  +
  scale_x_log10() +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)
df_full |>
  filter(country == "South Korea") |> 
  group_by(lat, long) |> 
  mutate(
    median_doy = median(bloom_doy, na.rm = TRUE)
  ) |> 
  ungroup() |> 
  distinct(lat, long, .keep_all = TRUE)  |> 
  ggplot(aes(x = alt, y = bio12, color = bloom_doy)) +
  geom_point() +
  scale_color_viridis_c(
    trans = "log10",
    breaks = trans_breaks(trans = "log10",
                          inv = function(x) round(10 ^ x, digits = 1))
  ) +
  labs(y = "Annual precipitation (mm)",
       x = "Altitude",
       color = "DOY (median)",
       title = "South Korea") +
  theme(legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm")) +
  geom_smooth(method = "gam",
              formula = y ~ ns(x, df = 2),
              color = "red",
              size = 0.5)

Interaction

En este caso como las variables son numéricas simplemento multiplico las dos variables bajo análisis para evaluar el efecto de interacción. Teniendo en cuenta que bloom_doy es la variable respuesta $(y)$ y si ejemplificamos un modelo sencillo que incluye las variables temperatura anual $(x_1)$ y precipitación $(x_2)$ como predictoras, mi intención es validar si el efecto acumulativo o conjunto está presente. Podemos estimar el efecto marginal de $x_1$ y definirlo como $\beta_1$, lo propio podemos hacer con el efecto marginal de $x_2$ y definirlo como $\beta_2$, sin embargo, también es de interés evaluar el efecto de interacción $(\beta_3)x_1x_2$, en ese orden de ideas se puede expresar el modelo matamático como se muestra a continuación:

\[y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_1x_2 + error\]

Un primer resultado interesante de las interaccione evaluadas es que el comportamiento discrepa entre países, por ejemplo, la interacción temperatura-precipitación muestra una relación inversa con bloom_doy, pero no es lineal en los tres países, en Corea del Sur la relación tiende a ser cuadrática. Algo similar podemos observar en la interacción temperatura-nitrógeno, en Japón y Corea del Sur no exhibe ninguna relación con la variable objetivo, sin embargo, en Suiza se evidencia una relación lineal inversa. En los tres países la interacción temperatura-bulk density se relaciona de forma negativa con la variable respuesta, pero en Japón es evidente la dependencia lineal que existe entre la variable $y$ y la predictora derivada de esta interacción.

Más adelante (en el documento 07-FeatureEngineering.qmd) hago uso de la regresión lasso para seleccionar predictores y evaluar posibles efectos de interacción que no fueron representados en este análisis exploratorio.

Code

countries <- c("Japan", "Switzerland", "South Korea")

df_full |>
    filter(country == countries[1]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * bio12) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Annual precipitation (mm)",
         y = "DOY",
         title = countries[1])
df_full |>
    filter(country == countries[2]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * bio12) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Annual precipitation (mm)",
         y = "DOY",
         title = countries[2])
df_full |>
    filter(country == countries[3]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * bio12) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Annual precipitation (mm)",
         y = "DOY",
         title = countries[3])

Code

df_full |>
    filter(country == countries[1]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `nitrogen_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Nitrogen (cg/kg)",
         y = "DOY",
         title = countries[1])
df_full |>
    filter(country == countries[2]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `nitrogen_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Nitrogen (cg/kg)",
         y = "DOY",
         title = countries[2])
df_full |>
    filter(country == countries[3]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `nitrogen_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Nitrogen (cg/kg)",
         y = "DOY",
         title = countries[3])

Code

df_full |>
    filter(country == countries[1]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `soc_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Soil organic carbon (dg/kg)",
         y = "DOY",
         title = countries[1])
df_full |>
    filter(country == countries[2]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `soc_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Soil organic carbon (dg/kg)",
         y = "DOY",
         title = countries[2])
df_full |>
    filter(country == countries[3]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `soc_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Soil organic carbon (dg/kg)",
         y = "DOY",
         title = countries[3])

Code

df_full |>
    filter(country == countries[1]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `bdod_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Bulk density (cg/cm³)",
         y = "DOY",
         title = countries[1])
df_full |>
    filter(country == countries[2]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `bdod_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Bulk density (cg/cm³)",
         y = "DOY",
         title = countries[2])
df_full |>
    filter(country == countries[3]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio1 * `bdod_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual mean temperature (°C) * Bulk density (cg/cm³)",
         y = "DOY",
         title = countries[3])

Code

df_full |>
    filter(country == countries[1]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio12 * `nitrogen_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual precipitation (mm) * Nitrogen (cg/kg)",
         y = "DOY",
         title = countries[1])
df_full |>
    filter(country == countries[2]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio12 * `nitrogen_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual precipitation (mm) * Nitrogen (cg/kg)",
         y = "DOY",
         title = countries[2])
df_full |>
    filter(country == countries[3]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = bio12 * `nitrogen_0-5cm_mean`) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Annual precipitation (mm) * Nitrogen (cg/kg)",
         y = "DOY",
         title = countries[3])

Code

df_full |>
    filter(country == countries[1]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = alt * bio12) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Altitude * Annual precipitation (mm) ",
         y = "DOY",
         title = countries[1])
df_full |>
    filter(country == countries[2]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = alt * bio12) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Altitude * Annual precipitation (mm) ",
         y = "DOY",
         title = countries[2])
df_full |>
    filter(country == countries[3]) |>
    group_by(lat, long) |>
    mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
    ungroup() |>
    distinct(lat, long, .keep_all = TRUE)  |>
    mutate(var_inter = alt * bio12) |>
    ggplot(aes(x = var_inter, y = bloom_doy)) +
    geom_point() +
    scale_x_log10() +
    geom_smooth(
      method = "gam",
      formula = y ~ ns(x, df = 2),
      color = "red",
      size = 0.5
    ) +
    labs(x = "Altitude * Annual precipitation (mm) ",
         y = "DOY",
         title = countries[3])

Principal Component Analysis (PCA)

En estos gráficos aplico la función distinct() por latitud y longitud porque la misma coordenada (así sean en años diferentes) tiene el mismo dato de bioclimáticas, suelo y coberturas, ya que esta información no posee temporalidad. Por esa razón resumo con la mediana del bloom_doy (ya que esta variable sí cambia con el tiempo).
Aprovechando que existe alta correlación entre las variables predictoras, ejecuto análisis de componentes principales para reducir la dimensionalidad e intentar ver si en un espacio más pequeño de variables se exhibe algún patrón de comportamiento relacionado con la variable respuesta.
Para facilitar las visualizaciones discretizo la variable respuesta en los siguientes niveles (ordinales):
- Q1: valores de bloom_doy inferiores al valor del cuartil 1: $DOY < Q1$.
- Q2: valores de bloom_doy mayores o iguales al valor del cuartil 1 y menores al valor del cuartil 2: $Q1 \leq DOY < Q2$.
- Q3: valores de bloom_doy mayores o iguales al valor del cuartil 2 y menores al valor del cuartil 3: $Q2 \leq DOY < Q3$.
- Q4: valores de bloom_doy superiores o iguales al valor del cuartil 3: $DOY \geq Q3$.
La variable respuesta discretizada doy_categ la introduzco al PCA como variable suplementaría cualitativa.
Independientemente del país las floraciones que no superan el Q1 (DOY < Q1) y las que se tardan más del Q3 (DOY >= Q3) exhiben diferencias al graficarlas en los tres primeros componentes principales, quedando ubicadas en lugares diferentes. Es importante mencionar que el patrón es menos visible en Corea del Sur. Floraciones que están entre el el Q1 y Q3 tienden a ser similares. Esta diferenciación que se aprecia entre floraciones precoces y floraciones tardías podría constituirse en un indicativo de perfiles diferenciales en el nicho ecológico o entorno ambiental al cual están expuestos los cerezos.
En el caso de Japón se observa que el componente 1 está asociado de forma positiva con variables que recogen información de temperatura (bio1, bio5, bio10, bio11) y precipitación (bio12, bio13, bio16), indicando que valores mayores a cero en el componente 1 los podemos asociar con valores altos de temperatura y precipitación, justamente bajo estas condiciones observamos las floraciones rápidas (puntos azules), de otro lado, este componente está asociado de forma negativa con la latitud, la longitud, bio4 (estacionalidad de la temperatura), bio7 (rango anual de temperatura) y en menor medida con el nitrógeno, indicando que valores mayores a cero de este componente estarán relacionados no sólo con valores altos de temperatura y precipitación sino también con menor variación en la temperatura (bio4) y menor rango de temperatura anual (bio7), por tanto, el perfil descrito es el que aplica a las floraciones precoces, justo lo contrario ocurre con floraciones tardías, donde se esperan temperaturas y precipitaciones bajas, con mayores variaciones en la temperatura a través del año, también podríamos afirmar que plantas que están ubicadas más al sur de Japón tienden a presentar floraciones retardadas. Por la ubicación de los puntos también podríamos afirmar que las floraciones precoces en Japón están moduladas por cambios de temperatura, sin embargo, otros factores como la capacidad de intercambio catiónico (cec) y la densidad del carbono orgánico (ocd) podrían ser influyentes, ya que se observa que estas variables tiene correlación positiva con el componente 2, indicando que valores mayores a cero en el componente 2 están asociados con suelos altos en capacidad de intercambio catiónico y suelos con mayor densidad del carbono orgánico.
En Suiza también se evidencia la diferencia entre floraciones precoces y floraciones tardías, no obstante, la correlación de los componentes con las variables permite inferir patrones de comportamiento ligeramente diferentes a Japón y Corea del Sur. El componente uno tiene correlación positiva con variables que recogen información de temperatura (bio4, bio5, bio10, bio11) y exhibe correlación negativa con la altitud, el contenido de carbono orgánico del suelo (soc), el inventario de carbono orgánico del suelo (ocs) y en menor medida con variables de precipitación (bio12, bio16, bio17), estas correlaciones permiten inferior que floraciones precoses se dan en condiciones de alta temperatura (también la variación anual) y en sitios donde hay menor cantidad de lluvia anual, menor contenido de carbono orgánico del suelo y ubicadas en altitudes bajas. Los cerezos que exhiben floraciones tardías tienden a estar ubicados en altitudes elevadas, sitios con mayor cantidad de precipitación anual, mayores contenidos de carbono orgánico del suelo.
En Corea del Sur se observa un comportamiento similar al de Japón y Suiza en cuanto a que la temperatura es el factor determinante en los tiempos de floración. Sin embargo, el componente 1 está asociado de forma negativa con variables que recogen información de temperatura (bio1, bio6, bio9, bio11) y de forma positiva con variables que recogen información de la variación anual de la temperatura (bio4, bio7). Las floraciones tardías se ubican a la derecha del 0 en el componente 1, es decir, que estas plantas se espera estén en lugares con mayores variaciones en la temperatura a través del año y también se espera que las temperaturas de estos sitios sean más bajas. Basado en la asociación positiva que tiene el componente 3 con las precipitaciones (bio12, bio15, bio16) y observando la ubicación de los puntos verdes en el gráfico de CP1 vs CP3, es posible intuir que las floraciones demoradas también están asociadas con mayores precipitaciones.
Es importante mencionar que todas estas interpretaciones las hago desde la perspectiva exploratoria, siguiendo las premisas de J. W. Tukey como se muestra en el siguiente diagrama. La idea principal ha sido reconocer patrones que me permitan realizar un diagnóstico objetivo de cara a la fase de modelación.

Con cuatro componentes principales se retiene el 69.4% de la variabilidad total.
Con diez componentes principales se retiene el 89.2% de la variabilidad total.
Disminuir de 41 dimensiones a 4 perdiente aproximadamente el 30% de la variabilidad no está nada mal para propósitos de graficación, sin embargo, para propósitos de modelación, valdría la pena reducir de 41 variables a solo 10 componentes.

Code

df_pca_japan1 <-
  df_full |>
  filter(country == countries[1]) 

value_q1 <- 
  quantile(df_pca_japan1$bloom_doy, probs = 0.25, na.rm = TRUE)

value_q2 <- 
  quantile(df_pca_japan1$bloom_doy, probs = 0.50, na.rm = TRUE)

value_q3 <- 
  quantile(df_pca_japan1$bloom_doy, probs = 0.75, na.rm = TRUE)

order_categ <-
  c("DOY < Q1",
    "Q1 <= DOY < Q2",
    "Q2 <= DOY < Q3",
    "DOY >= Q3")

df_pca_japan2 <-
  df_pca_japan1 |>
  group_by(lat, long) |>
  mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
  ungroup() |>
  distinct(lat, long, .keep_all = TRUE) |>
  select(-shrubs) |>
  mutate(
    doy_categ = case_when(
      median_doy < value_q1 ~ "DOY < Q1",
      median_doy >= value_q1 &
        median_doy < value_q2 ~ "Q1 <= DOY < Q2",
      median_doy >= value_q2 &
        median_doy < value_q3 ~ "Q2 <= DOY < Q3",
      median_doy >= value_q3 ~ "DOY >= Q3",
    ),
    doy_categ = factor(doy_categ, levels = order_categ)
  ) |>
  select(-c(location, country, bloom_date, bloom_doy, median_doy)) |>
  relocate(doy_categ, everything())


pca_japan <- PCA(X = df_pca_japan2, graph = FALSE, quali.sup = 1)

df_pca_japan2$pc1 <- pca_japan$ind$coord[, 1]
df_pca_japan2$pc2 <- pca_japan$ind$coord[, 2]
df_pca_japan2$pc3 <- pca_japan$ind$coord[, 3]

eigen_pc1 <- pca_japan$eig[, 2][1] |> round(digits = 1)
eigen_pc2 <- pca_japan$eig[, 2][2] |> round(digits = 1)
eigen_pc3 <- pca_japan$eig[, 2][3] |> round(digits = 1)

df_pca_japan2 |> 
  ggplot(aes(x = pc1, y = pc2, color = doy_categ)) +
  geom_point(size = 2.5, alpha = 0.75, shape = 18) +
  geom_hline(yintercept = 0, lty = 2, color = "firebrick3") +
  geom_vline(xintercept = 0, lty = 2, color = "firebrick3") +
  labs(x = glue("Component 1 ({eigen_pc1}%)"),
       y = glue("Component 2 ({eigen_pc2}%)")) +
  scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) +
  labs(color = "", title = "Component 1 vs Component 2")
df_pca_japan2 |>
  select(-c(doy_categ)) |> 
  correlate(method = "pearson") |> 
  filter(term %in% c("pc1", "pc2", "pc3")) |> 
  select(-c(pc1, pc2, pc3)) |> 
  pivot_longer(-term) |> 
  mutate(
    name = str_replace_all(
      name,
      "wildareas-v3-1993-human-footprint",
      "human-footprint93"
    ),
    name = str_replace_all(
      name,
      "wildareas-v3-2009-human-footprint",
      "human-footprint09"
    )
  ) |> 
  ggplot(aes(x = reorder(name, value), y = term, fill = value)) +
  geom_tile() + 
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  labs(x = "", y = "Component", fill = "Correlation") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7.5),
        legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))
df_pca_japan2 |> 
  ggplot(aes(x = pc1, y = pc3, color = doy_categ)) +
  geom_point(size = 2.5, alpha = 0.75, shape = 18) +
  geom_hline(yintercept = 0, lty = 2, color = "firebrick3") +
  geom_vline(xintercept = 0, lty = 2, color = "firebrick3") +
  labs(x = glue("Component 1 ({eigen_pc1}%)"),
       y = glue("Component 3 ({eigen_pc3}%)")) +
  scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) +
  labs(color = "", title = "Component 1 vs Component 3")

Con cuatro componentes principales se retiene el 66.4% de la variabilidad total.
Con diez componentes principales se retiene el 85.5% de la variabilidad total.

Code

df_pca_switzerland1 <-
  df_full |>
  filter(country == countries[2]) 

value_q1 <- 
  quantile(df_pca_switzerland1$bloom_doy, probs = 0.25, na.rm = TRUE)

value_q2 <- 
  quantile(df_pca_switzerland1$bloom_doy, probs = 0.50, na.rm = TRUE)

value_q3 <- 
  quantile(df_pca_switzerland1$bloom_doy, probs = 0.75, na.rm = TRUE)

order_categ <-
  c("DOY < Q1",
    "Q1 <= DOY < Q2",
    "Q2 <= DOY < Q3",
    "DOY >= Q3")

df_pca_switzerland2 <-
  df_pca_switzerland1 |>
  group_by(lat, long) |>
  mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
  ungroup() |>
  distinct(lat, long, .keep_all = TRUE) |>
  select(-shrubs) |>
  mutate(
    doy_categ = case_when(
      median_doy < value_q1 ~ "DOY < Q1",
      median_doy >= value_q1 &
        median_doy < value_q2 ~ "Q1 <= DOY < Q2",
      median_doy >= value_q2 &
        median_doy < value_q3 ~ "Q2 <= DOY < Q3",
      median_doy >= value_q3 ~ "DOY >= Q3",
    ),
    doy_categ = factor(doy_categ, levels = order_categ)
  ) |>
  select(-c(location, country, bloom_date, bloom_doy, median_doy)) |>
  relocate(doy_categ, everything())


pca_switzerland <- PCA(X = df_pca_switzerland2, graph = FALSE, quali.sup = 1)

df_pca_switzerland2$pc1 <- pca_switzerland$ind$coord[, 1]
df_pca_switzerland2$pc2 <- pca_switzerland$ind$coord[, 2]
df_pca_switzerland2$pc3 <- pca_switzerland$ind$coord[, 3]

eigen_pc1 <- pca_switzerland$eig[, 2][1] |> round(digits = 1)
eigen_pc2 <- pca_switzerland$eig[, 2][2] |> round(digits = 1)
eigen_pc3 <- pca_switzerland$eig[, 2][3] |> round(digits = 1)

df_pca_switzerland2 |> 
  ggplot(aes(x = pc1, y = pc2, color = doy_categ)) +
  geom_point(size = 2.5, alpha = 0.75, shape = 18) +
  geom_hline(yintercept = 0, lty = 2, color = "firebrick3") +
  geom_vline(xintercept = 0, lty = 2, color = "firebrick3") +
  labs(x = glue("Component 1 ({eigen_pc1}%)"),
       y = glue("Component 2 ({eigen_pc2}%)")) +
  scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) +
  labs(color = "", title = "Component 1 vs Component 2")
df_pca_switzerland2 |>
  select(-c(doy_categ)) |> 
  correlate(method = "pearson") |> 
  filter(term %in% c("pc1", "pc2", "pc3")) |> 
  select(-c(pc1, pc2, pc3)) |> 
  pivot_longer(-term) |> 
  mutate(
    name = str_replace_all(
      name,
      "wildareas-v3-1993-human-footprint",
      "human-footprint93"
    ),
    name = str_replace_all(
      name,
      "wildareas-v3-2009-human-footprint",
      "human-footprint09"
    )
  ) |> 
  ggplot(aes(x = reorder(name, value), y = term, fill = value)) +
  geom_tile() + 
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  labs(x = "", y = "Component", fill = "Correlation") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7.5),
        legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))
df_pca_switzerland2 |> 
  ggplot(aes(x = pc1, y = pc3, color = doy_categ)) +
  geom_point(size = 2.5, alpha = 0.75, shape = 18) +
  geom_hline(yintercept = 0, lty = 2, color = "firebrick3") +
  geom_vline(xintercept = 0, lty = 2, color = "firebrick3") +
  labs(x = glue("Component 1 ({eigen_pc1}%)"),
       y = glue("Component 3 ({eigen_pc3}%)")) +
  scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) +
  labs(color = "", title = "Component 1 vs Component 3")

Con cuatro componentes principales se retiene el 66.4% de la variabilidad total.
Con diez componentes principales se retiene el 89.7% de la variabilidad total.

Code

df_pca_southk1 <-
  df_full |>
  filter(country == countries[3]) 

value_q1 <- 
  quantile(df_pca_southk1$bloom_doy, probs = 0.25, na.rm = TRUE)

value_q2 <- 
  quantile(df_pca_southk1$bloom_doy, probs = 0.50, na.rm = TRUE)

value_q3 <- 
  quantile(df_pca_southk1$bloom_doy, probs = 0.75, na.rm = TRUE)

order_categ <-
  c("DOY < Q1",
    "Q1 <= DOY < Q2",
    "Q2 <= DOY < Q3",
    "DOY >= Q3")

df_pca_southk2 <-
  df_pca_southk1 |>
  group_by(lat, long) |>
  mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |>
  ungroup() |>
  distinct(lat, long, .keep_all = TRUE) |>
  select(-shrubs) |>
  mutate(
    doy_categ = case_when(
      median_doy < value_q1 ~ "DOY < Q1",
      median_doy >= value_q1 &
        median_doy < value_q2 ~ "Q1 <= DOY < Q2",
      median_doy >= value_q2 &
        median_doy < value_q3 ~ "Q2 <= DOY < Q3",
      median_doy >= value_q3 ~ "DOY >= Q3",
    ),
    doy_categ = factor(doy_categ, levels = order_categ)
  ) |>
  select(-c(location, country, bloom_date, bloom_doy, median_doy)) |>
  relocate(doy_categ, everything())


pca_southk <- PCA(X = df_pca_southk2, graph = FALSE, quali.sup = 1)

df_pca_southk2$pc1 <- pca_southk$ind$coord[, 1]
df_pca_southk2$pc2 <- pca_southk$ind$coord[, 2]
df_pca_southk2$pc3 <- pca_southk$ind$coord[, 3]

eigen_pc1 <- pca_southk$eig[, 2][1] |> round(digits = 1)
eigen_pc2 <- pca_southk$eig[, 2][2] |> round(digits = 1)
eigen_pc3 <- pca_southk$eig[, 2][3] |> round(digits = 1)

df_pca_southk2 |> 
  ggplot(aes(x = pc1, y = pc2, color = doy_categ)) +
  geom_point(size = 2.5, alpha = 0.75, shape = 18) +
  geom_hline(yintercept = 0, lty = 2, color = "firebrick3") +
  geom_vline(xintercept = 0, lty = 2, color = "firebrick3") +
  labs(x = glue("Component 1 ({eigen_pc1}%)"),
       y = glue("Component 2 ({eigen_pc2}%)")) +
  scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) +
  labs(color = "", title = "Component 1 vs Component 2")
df_pca_southk2 |>
  select(-c(doy_categ)) |> 
  correlate(method = "pearson") |> 
  filter(term %in% c("pc1", "pc2", "pc3")) |> 
  select(-c(pc1, pc2, pc3)) |> 
  pivot_longer(-term) |> 
  mutate(
    name = str_replace_all(
      name,
      "wildareas-v3-1993-human-footprint",
      "human-footprint93"
    ),
    name = str_replace_all(
      name,
      "wildareas-v3-2009-human-footprint",
      "human-footprint09"
    )
  ) |> 
  ggplot(aes(x = reorder(name, value), y = term, fill = value)) +
  geom_tile() + 
  scale_fill_gradient2(
    low = muted("red"),
    mid = "white",
    high = muted("dodgerblue2"),
    midpoint = 0
  ) +
  labs(x = "", y = "Component", fill = "Correlation") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7.5),
        legend.key.size = unit(1, "cm"),
        legend.key.width = unit(2, "cm"))
df_pca_southk2 |> 
  ggplot(aes(x = pc1, y = pc3, color = doy_categ)) +
  geom_point(size = 2.5, alpha = 0.75, shape = 18) +
  geom_hline(yintercept = 0, lty = 2, color = "firebrick3") +
  geom_vline(xintercept = 0, lty = 2, color = "firebrick3") +
  labs(x = glue("Component 1 ({eigen_pc1}%)"),
       y = glue("Component 3 ({eigen_pc3}%)")) +
  scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) +
  labs(color = "", title = "Component 1 vs Component 3")

--- title: "International Cherry Blossom Prediction Competition" subtitle: "Exploratory Data Analysis 1" author: "Edimer David Jaramillo" date: "`r lubridate::now()`" lang: en-US format: html: page-layout: article toc: true code-fold: true df-print: paged toc-location: left number-depth: 4 theme: yeti code-copy: true highlight-style: github embed-resources: true code-tools: source: true --- ```{r} #| label: setup #| include: false knitr::opts_chunk$set(echo = TRUE, warning = FALSE, error = FALSE, message = FALSE, fig.align = 'center') ``` # Document description - En este documento se lleva a cabo análisis exploratorio de datos con las variables bioclimáticas, de suelo, coberturas terrestres y huella humana. Recordemos que las variables bioclimáticas son un consolidado de los años 1970 a 2000. Este documento tiene como finalidad comparar las fechas de floración entre países, validar relaciones que han sido reportadas en literatura (por ejemplo, la relación de temperatura con *DOY*), determinar la existencia de efectos de interacción y establecer cuáles de las variables predictoras podrían ser eventualmente candidatas para incluir en la modelación (*selección de predictores*). Podríamos decir que las variables analizadas en este documento son *"estáticas"* o *"fijas"*, ya que no tenemos temporalidad de ellas, son datos consolidados en ventanas de tiempo preestablecidas. - En el siguiente documento (**06-EDA2.qmd**) se muestran resultados de análisis exploratorio para variables *"dinámicas"*, me refiero a ellas como dinámicas ya que las variables climatológicas y de fotoperíodo tienen temporalidad (datos diarios) y obtuve información desde el 01 de enero de 1981 hasta el 31 de diciembre de 2023. # Libraries and setup ```{r} library(tidyverse) library(arrow) library(splines) library(corrr) library(scales) library(FactoMineR) library(factoextra) library(glue) theme_set(theme_bw() + theme(legend.position = "top")) ``` # Data - Extraigo una variable de nombre `country` que tiene el país, esto lo hago con la finalidad de evaluar patrones de comportamiento por diferentes países y porque podría ser de ayuda para extrapolar a regiones donde no tenemos registros históricos de floración. Para diferenciar los datos de *NPN* para USA, los etiqueto en la variable `country` como `USA-NPN` y a los de Washington DC como `USA-WDC`. - En varios de los resultados de ese documento no se incluye `USA-WDC` porque sólo es una coordenada (lat=38.88535, long=-77.03863) y tampoco se muestran los resultados de `USA-NPN`, la razón es que las variables bioclimáticas no están para las coordenadas de `USA-NPN`, sin embargo, como información de clima fue suministrada por *NPN*, uso esta información para análisis exploratorio al final de este documento y también en el próximo (**06-EDA2.qmd**). ```{r} df_full <- read_parquet("../external-data/df_full.parquet") |> mutate( country = case_when( str_detect(location, "Japan") ~ "Japan", str_detect(location, "kyoto") ~ "Japan", str_detect(location, "liestal") ~ "Switzerland", str_detect(location, "Switzerland") ~ "Switzerland", str_detect(location, "South Korea") ~ "South Korea", str_detect(location, "vancouver") ~ "Canada", str_detect(location, "washingtondc") ~ "USA-WDC", str_detect(location, "NPN") ~ "USA-NPN" ) ) |> relocate(location, country, everything()) ``` # DOY distribution by country ```{r} #| fig-width: 7 #| fig-height: 5.5 #| layout-nrow: 1 #| column: screen-inset-shaded df_full |> filter(location != "NPN") |> ggplot(aes(x = bloom_doy)) + facet_wrap(~country, ncol = 1, scales = "free_y") + geom_histogram(color = "black") + geom_rug() + labs(x = "Bloom DOY", y = "Count", title = "DOY distribution by country", subtitle = "Original scale") df_full |> filter(location != "NPN") |> ggplot(aes(x = bloom_doy)) + facet_wrap(~country, ncol = 1, scales = "free_y") + geom_histogram(color = "black") + geom_rug() + scale_x_log10() + labs(x = "Bloom DOY", y = "Count", title = "DOY distribution by country", subtitle = "Logarithmic scale") ``` # Scatter plot of all vs DOY - Los siguientes gráficos permiten evidenciar relaciones (lineales y no lineales) entre algunas de las variables bioclimáticas y el día de la floración, también se observan asociaciones entre composición de suelo, por ejemplo, pH y *DOY*, no son dicientes o claras las relaciones entre las variables de cobertura terreste, huella humana, pastizales, arbutos, agua y la fecha de la floración. - Con estos gráficos se valida o confirma que la temperatura (*bio1*, *bio5*, *bio8*, *bio9*, *bio10* y *bio11*) y el día de la floración exhiben una relación que podríamos sugerir está entre lineal y cuadrática, independiente del país, a mayor temperatura se han observado floraciones más precoces, este comportamiento es similar en Japón y Suiza, no obstante, Corea del Sur presenta algunas diferencias, por ejemplo, en la variable *bio11* para Japón y Suiza, la relación de tipo cuadrático es evidente, pero en Corea del sur no lo es. - Las precipitaciones (*bio12* a *bio19*) también muestran relaciones interesantes con la fecha de floración, con ligeras diferencias entre países, por ejemplo, la precipitación anual (*bio12*) muestra una relación lineal inversa con el *DOY*, en Suiza por el contrario observamos una relación cuadrática positiva y en Corea del Sur no es clara la relación. Este resultado puede sugerir que temperatura y precipitación son factores que interactuan y tener un efecto conjunto; más adelante se exploran interacciones entre variables para determinar covariaciones de interés. - La estacionalidad de la temperatura (*bio4*) muestra un patrón de comportamiento diferente entre países, en Japón la relación es positiva y lineal, en Corea del Sur no es clara la asociación y en Suiza muesra una relación inversa cuadrática. Recordemos que las variables bioclimáticas (*bio_*) recogen información de los años 1970 a 2000, con una resolución de $1 km^2$, este hallazgo será comparado más adelante con la información diaria de temperatura para validar las relaciones que se evidencian en este gráfico. - En cuanto a la composición del suelo, el nitrógeno, el pH del agua del suelo, el contenido de carbono orgánico del suelo (*SOC*), la densidad aparente (*BDOD*) y la capacidad de intercambio catiónico (*CEC*) muestran tendencias en la relación con el *DOY*. Las otras variables parecen no tener ninguna asociación (lineal o no lineal) con la variable objetivo. - Finalmente, la información de tierras de cultivo (*CROPLAND*), huella humana (*FOOTPRINT*) y coberturas terrestres, no exhiben ningún tipo de relación que pudiera ser considerada de interés para análisis posteriores, además, algunas de estas variables (por ejemplo *SHRUBS*) carecen de información para las coordenadas bajo análisis. ::: {.panel-tabset} ## Japan ```{r} #| fig-width: 9 #| fig-height: 16 #| fig-align: center #| column: page-inset-right df_full |> filter(country == "Japan") |> select(-c(location, country, year, bloom_date)) |> pivot_longer(cols = -bloom_doy) |> ggplot(aes(x = value, y = bloom_doy)) + facet_wrap( ~ name, scales = "free", ncol = 4) + geom_point(size = 0.5, alpha = 0.25) + stat_bin_2d(aes(fill = ..density..), geom = "raster", contour = FALSE, show.legend = FALSE) + scale_fill_distiller(palette = 4, direction = -1) + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "firebrick3", se = FALSE) + scale_x_log10() + labs(y = "DOY", x = "", title = "Japan") ``` ## Switzerland ```{r} #| fig-width: 9 #| fig-height: 16 #| fig-align: center #| column: page-inset-right df_full |> filter(country == "Switzerland") |> select(-c(location, country, year, bloom_date)) |> pivot_longer(cols = -bloom_doy) |> ggplot(aes(x = value, y = bloom_doy)) + facet_wrap( ~ name, scales = "free", ncol = 4) + geom_point(size = 0.5, alpha = 0.25) + stat_bin_2d(aes(fill = ..density..), geom = "raster", contour = FALSE, show.legend = FALSE) + scale_fill_distiller(palette = 4, direction = -1) + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "firebrick3", se = FALSE) + scale_x_log10() + labs(y = "DOY", x = "", title = "Switzerland") ``` ## South Korea ```{r} #| fig-width: 9 #| fig-height: 16 #| fig-align: center #| column: page-inset-right df_full |> filter(country == "South Korea") |> select(-c(location, country, year, bloom_date)) |> pivot_longer(cols = -bloom_doy) |> ggplot(aes(x = value, y = bloom_doy)) + facet_wrap( ~ name, scales = "free", ncol = 4) + geom_point(size = 0.5, alpha = 0.25) + stat_bin_2d(aes(fill = ..density..), geom = "raster", contour = FALSE, show.legend = FALSE) + scale_fill_distiller(palette = 4, direction = -1) + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "firebrick3", se = FALSE) + scale_x_log10() + labs(y = "DOY", x = "", title = "South Korea") ``` ## All the countries ```{r} #| fig-width: 9 #| fig-height: 16 #| fig-align: center #| column: page-inset-right df_full |> filter(location != "NPN") |> select(-c(location, country, year, bloom_date)) |> pivot_longer(cols = -bloom_doy) |> ggplot(aes(x = value, y = bloom_doy)) + facet_wrap( ~ name, scales = "free", ncol = 4) + geom_point(size = 0.5, alpha = 0.25) + stat_bin_2d(aes(fill = ..density..), geom = "raster", contour = FALSE, show.legend = FALSE) + scale_fill_distiller(palette = 4, direction = -1) + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "firebrick3", se = FALSE) + scale_x_log10() + labs(y = "DOY", x = "", title = "All the countries") ``` ::: # Non-parametric correlations - Estos perfiles de correlación muestran para cada país y todos los países, de mayor a menor correlación (Spearman). En Japón y Corea del Sur, los cerezos ubicados más al norte parecen tardarse más en la floración (colores azules), la variable *latitud* es la primera en el perfil de correlaciones de estos países, sin embargo, en Suiza no es la latitud el factor lineal de mayor influcencia, es la *altitud* seguido de variables de composición de suelo y en el top 10 están tres variables que recogen información de precipitación. Si observamos la parte baja del gráfico (colores rojos) en todos los países coincide la variable temperatura anual (*bio1*) en ser el factor lineal inverso de mayor magnitud, seguido de otras variables indicadoras también de temperaturas. Este gráfico es de mucha utilidad no sólo para representar relaciones de tipo lineal sino también para obtener una caracterización rápida de la relación del clima con la fenología en cerezos ubicados en diferentes países, podemos decir que si pensamos en los factores que se asociación con **tardanza** en la floración, los perfiles entre países son muy diferentes, no obstante, si pensamos en los factores que se asociación con la **precocidad** en la floración, los perfiles entre países son muy parecidos. En conclusión, florecer más rápido parece tener unas *"causas"* comunes y globales ligadas a la temperatura, sin embargo, el retraso en la floración podría ser más difícil de establecer y generalizar. ```{r} #| fig-width: 4 #| fig-height: 8 #| column: screen #| layout-nrow: 1 df_full |> filter(country == "Japan") |> select(where(is.numeric), -c(shrubs)) |> correlate(method = "spearman") |> filter(term == "bloom_doy") |> pivot_longer(cols = -term) |> filter(!is.na(value)) |> mutate( name = str_replace_all( name, "wildareas-v3-1993-human-footprint", "human-footprint93" ), name = str_replace_all( name, "wildareas-v3-2009-human-footprint", "human-footprint09" ) ) |> ggplot(aes( x = reorder(name, value), y = term, fill = value )) + geom_tile() + scale_fill_gradient2( low = muted("red"), mid = "white", high = muted("dodgerblue2"), midpoint = 0 ) + theme(axis.text.x = element_text(angle = 35, hjust = 1)) + labs(x = "", y = "", fill = "", title = "Japan") + coord_flip() df_full |> filter(country == "Switzerland") |> select(where(is.numeric), -c(shrubs)) |> correlate(method = "spearman") |> filter(term == "bloom_doy") |> pivot_longer(cols = -term) |> filter(!is.na(value)) |> mutate(name = str_replace_all(name, "wildareas-v3-1993-human-footprint", "human-footprint93"), name = str_replace_all(name, "wildareas-v3-2009-human-footprint", "human-footprint09")) |> ggplot(aes(x = reorder(name, value), y = term, fill = value)) + geom_tile() + geom_tile() + scale_fill_gradient2( low = muted("red"), mid = "white", high = muted("dodgerblue2"), midpoint = 0 ) + theme(axis.text.x = element_text(angle = 35, hjust = 1)) + labs(x = "", y = "", fill = "", title = "Switzerland") + coord_flip() df_full |> filter(country == "South Korea") |> select(where(is.numeric), -c(shrubs)) |> correlate(method = "spearman") |> filter(term == "bloom_doy") |> pivot_longer(cols = -term) |> filter(!is.na(value)) |> mutate(name = str_replace_all(name, "wildareas-v3-1993-human-footprint", "human-footprint93"), name = str_replace_all(name, "wildareas-v3-2009-human-footprint", "human-footprint09")) |> ggplot(aes(x = reorder(name, value), y = term, fill = value)) + geom_tile() + geom_tile() + scale_fill_gradient2( low = muted("red"), mid = "white", high = muted("dodgerblue2"), midpoint = 0 ) + theme(axis.text.x = element_text(angle = 35, hjust = 1)) + labs(x = "", y = "", fill = "", title = "South Korea") + coord_flip() df_full |> filter(location != "NPN") |> select(where(is.numeric), -c(shrubs)) |> correlate(method = "spearman") |> filter(term == "bloom_doy") |> pivot_longer(cols = -term) |> filter(!is.na(value)) |> mutate(name = str_replace_all(name, "wildareas-v3-1993-human-footprint", "human-footprint93"), name = str_replace_all(name, "wildareas-v3-2009-human-footprint", "human-footprint09")) |> ggplot(aes(x = reorder(name, value), y = term, fill = value)) + geom_tile() + geom_tile() + scale_fill_gradient2( low = muted("red"), mid = "white", high = muted("dodgerblue2"), midpoint = 0 ) + theme(axis.text.x = element_text(angle = 35, hjust = 1), legend.key.size = unit(0.85, "cm")) + labs(x = "", y = "", fill = "", title = "All the countries") + coord_flip() ``` # Interaction of variables - Con los siguientes gráficos pretendo reconocer si existen comportamientos de interacción, es decir, intento confirmar a partir de representaciones visuales si hay efectos conjuntos entre las variables predictoras. En vista de que el conjunto de variables predictoras es "grande" $(P = 41)$ generar todas las interacciones (dobles, triples o de mayor jerarquía) podría ser dispendioso y demorado. Las siguientes son las dos aproximaciones que implemento: - **1.** Inicialmente voy a graficar diagramas de dispersión con las variables que mostraron tendencias en los gráficos de dispersión y los perfiles de correlación; el color de los puntos está determinado por la variable respuesta `bloom_doy`. - **2.** Grafico la interacción doble para algunas variables de interés vs la variable respuesta. ## Scatter plot - En estos gráficos aplico la función distinct() por latitud y longitud porque la misma coordenada (así sean en años diferentes) tiene el mismo dato de bioclimáticas, suelo y coberturas, ya que esta información no posee temporalidad. Por esa razón resumo con la **mediana del `bloom_doy`** (ya que esta variable sí cambia con el tiempo). ::: {.panel-tabset} ### Temperature - Precipitation ```{r} #| fig-width: 6 #| fig-height: 4.5 #| column: screen #| layout-nrow: 1 df_full |> filter(country == "Japan") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio1, y = bio12, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual mean temperature (°C)", y = "Annual precipitation (mm)", color = "DOY (median)", title = "Japan") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) df_full |> filter(country == "Switzerland") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio1, y = bio12, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual mean temperature (°C)", y = "Annual precipitation (mm)", color = "DOY (median)", title = "Switzerland") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) df_full |> filter(country == "South Korea") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio1, y = bio12, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual mean temperature (°C)", y = "Annual precipitation (mm)", color = "DOY (median)", title = "South Korea") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) ``` ### Temperature - Nitrogen ```{r} #| fig-width: 6 #| fig-height: 4.5 #| column: screen #| layout-nrow: 1 df_full |> filter(country == "Japan") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio1, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual mean temperature (°C)", y = "Nitrogen (cg/kg)", color = "DOY (median)", title = "Japan") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + scale_x_log10() + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) df_full |> filter(country == "Switzerland") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio1, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual mean temperature (°C)", y = "Nitrogen (cg/kg)", color = "DOY (median)", title = "Switzerland") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + scale_x_log10() + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) df_full |> filter(country == "South Korea") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio1, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual mean temperature (°C)", y = "Nitrogen (cg/kg)", color = "DOY (median)", title = "South Korea") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) ``` ### Temperature - SOC ```{r} #| fig-width: 6 #| fig-height: 4.5 #| column: screen #| layout-nrow: 1 df_full |> filter(`soc_0-5cm_mean` > 0) |> filter(country == "Japan") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio1, y = `soc_0-5cm_mean`, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual mean temperature (°C)", y = "Soil organic carbon (dg/kg)", color = "DOY (median)", title = "Japan") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + scale_x_log10() + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) df_full |> filter(`soc_0-5cm_mean` > 0) |> filter(country == "Switzerland") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio1, y = `soc_0-5cm_mean`, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual mean temperature (°C)", y = "Soil organic carbon (dg/kg)", color = "DOY (median)", title = "Switzerland") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + scale_x_log10() + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) df_full |> filter(`soc_0-5cm_mean` > 0) |> filter(country == "South Korea") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio1, y = `soc_0-5cm_mean`, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual mean temperature (°C)", y = "Soil organic carbon (dg/kg)", color = "DOY (median)", title = "South Korea") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) ``` ### Temperature - BDOD ```{r} #| fig-width: 6 #| fig-height: 4.5 #| column: screen #| layout-nrow: 1 df_full |> filter(`bdod_0-5cm_mean` > 0) |> filter(country == "Japan") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio1, y = `bdod_0-5cm_mean`, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual mean temperature (°C)", y = "Bulk density (cg/cm³)", color = "DOY (median)", title = "Japan") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + scale_x_log10() + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) df_full |> filter(`bdod_0-5cm_mean` > 0) |> filter(country == "Switzerland") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio1, y = `bdod_0-5cm_mean`, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual mean temperature (°C)", y = "Bulk density (cg/cm³)", color = "DOY (median)", title = "Switzerland") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + scale_x_log10() + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) df_full |> filter(`bdod_0-5cm_mean` > 0) |> filter(country == "South Korea") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio1, y = `bdod_0-5cm_mean`, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual mean temperature (°C)", y = "Bulk density (cg/cm³)", color = "DOY (median)", title = "South Korea") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) ``` ### Precipitation - Nitrogen ```{r} #| fig-width: 6 #| fig-height: 4.5 #| column: screen #| layout-nrow: 1 df_full |> filter(`nitrogen_0-5cm_mean` > 0) |> filter(country == "Japan") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio12, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual precipitation (mm)", y = "Nitrogen (cg/kg)", color = "DOY (median)", title = "Japan") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + scale_x_log10() + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) df_full |> filter(`nitrogen_0-5cm_mean` > 0) |> filter(country == "Switzerland") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio12, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual precipitation (mm)", y = "Nitrogen (cg/kg)", color = "DOY (median)", title = "Switzerland") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + scale_x_log10() + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) df_full |> filter(`nitrogen_0-5cm_mean` > 0) |> filter(country == "South Korea") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = bio12, y = `nitrogen_0-5cm_mean`, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(x = "Annual precipitation (mm)", y = "Nitrogen (cg/kg)", color = "DOY (median)", title = "South Korea") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) ``` ### Altitude - Precipitation ```{r} #| fig-width: 6 #| fig-height: 4.5 #| column: screen #| layout-nrow: 1 df_full |> filter(country == "Japan") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = alt, y = bio12, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(y = "Annual precipitation (mm)", x = "Altitude", color = "DOY (median)", title = "Japan") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + scale_x_log10() + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) df_full |> filter(country == "Switzerland") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = alt, y = bio12, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(y = "Annual precipitation (mm)", x = "Altitude)", color = "DOY (median)", title = "Switzerland") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + scale_x_log10() + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) df_full |> filter(country == "South Korea") |> group_by(lat, long) |> mutate( median_doy = median(bloom_doy, na.rm = TRUE) ) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> ggplot(aes(x = alt, y = bio12, color = bloom_doy)) + geom_point() + scale_color_viridis_c( trans = "log10", breaks = trans_breaks(trans = "log10", inv = function(x) round(10 ^ x, digits = 1)) ) + labs(y = "Annual precipitation (mm)", x = "Altitude", color = "DOY (median)", title = "South Korea") + theme(legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) + geom_smooth(method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5) ``` ::: ## Interaction En este caso como las variables son numéricas simplemento multiplico las dos variables bajo análisis para evaluar el efecto de interacción. Teniendo en cuenta que *bloom_doy* es la variable respuesta $(y)$ y si ejemplificamos un modelo sencillo que incluye las variables temperatura anual $(x_1)$ y precipitación $(x_2)$ como predictoras, mi intención es validar si el efecto acumulativo o conjunto está presente. Podemos estimar el efecto marginal de $x_1$ y definirlo como $\beta_1$, lo propio podemos hacer con el efecto marginal de $x_2$ y definirlo como $\beta_2$, sin embargo, también es de interés evaluar el efecto de interacción $(\beta_3)x_1x_2$, en ese orden de ideas se puede expresar el modelo matamático como se muestra a continuación: $$y = \beta_0 + \beta_1x_1 + \beta_2x_2 + \beta_3x_1x_2 + error$$ Un primer resultado interesante de las interaccione evaluadas es que el comportamiento discrepa entre países, por ejemplo, la interacción *temperatura-precipitación* muestra una relación inversa con *bloom_doy*, pero no es lineal en los tres países, en Corea del Sur la relación tiende a ser cuadrática. Algo similar podemos observar en la interacción *temperatura-nitrógeno*, en Japón y Corea del Sur no exhibe ninguna relación con la variable objetivo, sin embargo, en Suiza se evidencia una relación lineal inversa. En los tres países la interacción *temperatura-bulk density* se relaciona de forma negativa con la variable respuesta, pero en Japón es evidente la dependencia lineal que existe entre la variable $y$ y la predictora derivada de esta interacción. Más adelante (en el documento **07-FeatureEngineering.qmd**) hago uso de la *regresión lasso* para seleccionar predictores y evaluar posibles efectos de interacción que no fueron representados en este análisis exploratorio. ::: {.panel-tabset} ### Temperature - Precipitation ```{r} #| fig-width: 6 #| fig-height: 4.5 #| column: screen #| layout-nrow: 1 countries <- c("Japan", "Switzerland", "South Korea") df_full |> filter(country == countries[1]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio1 * bio12) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual mean temperature (°C) * Annual precipitation (mm)", y = "DOY", title = countries[1]) df_full |> filter(country == countries[2]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio1 * bio12) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual mean temperature (°C) * Annual precipitation (mm)", y = "DOY", title = countries[2]) df_full |> filter(country == countries[3]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio1 * bio12) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual mean temperature (°C) * Annual precipitation (mm)", y = "DOY", title = countries[3]) ``` ### Temperature - Nitrogen ```{r} #| fig-width: 6 #| fig-height: 4.5 #| column: screen #| layout-nrow: 1 df_full |> filter(country == countries[1]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio1 * `nitrogen_0-5cm_mean`) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual mean temperature (°C) * Nitrogen (cg/kg)", y = "DOY", title = countries[1]) df_full |> filter(country == countries[2]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio1 * `nitrogen_0-5cm_mean`) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual mean temperature (°C) * Nitrogen (cg/kg)", y = "DOY", title = countries[2]) df_full |> filter(country == countries[3]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio1 * `nitrogen_0-5cm_mean`) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual mean temperature (°C) * Nitrogen (cg/kg)", y = "DOY", title = countries[3]) ``` ### Temperature - SOC ```{r} #| fig-width: 6 #| fig-height: 4.5 #| column: screen #| layout-nrow: 1 df_full |> filter(country == countries[1]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio1 * `soc_0-5cm_mean`) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual mean temperature (°C) * Soil organic carbon (dg/kg)", y = "DOY", title = countries[1]) df_full |> filter(country == countries[2]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio1 * `soc_0-5cm_mean`) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual mean temperature (°C) * Soil organic carbon (dg/kg)", y = "DOY", title = countries[2]) df_full |> filter(country == countries[3]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio1 * `soc_0-5cm_mean`) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual mean temperature (°C) * Soil organic carbon (dg/kg)", y = "DOY", title = countries[3]) ``` ### Temperature - BDOD ```{r} #| fig-width: 6 #| fig-height: 4.5 #| column: screen #| layout-nrow: 1 df_full |> filter(country == countries[1]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio1 * `bdod_0-5cm_mean`) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual mean temperature (°C) * Bulk density (cg/cm³)", y = "DOY", title = countries[1]) df_full |> filter(country == countries[2]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio1 * `bdod_0-5cm_mean`) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual mean temperature (°C) * Bulk density (cg/cm³)", y = "DOY", title = countries[2]) df_full |> filter(country == countries[3]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio1 * `bdod_0-5cm_mean`) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual mean temperature (°C) * Bulk density (cg/cm³)", y = "DOY", title = countries[3]) ``` ### Precipitation - Nitrogen ```{r} #| fig-width: 6 #| fig-height: 4.5 #| column: screen #| layout-nrow: 1 df_full |> filter(country == countries[1]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio12 * `nitrogen_0-5cm_mean`) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual precipitation (mm) * Nitrogen (cg/kg)", y = "DOY", title = countries[1]) df_full |> filter(country == countries[2]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio12 * `nitrogen_0-5cm_mean`) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual precipitation (mm) * Nitrogen (cg/kg)", y = "DOY", title = countries[2]) df_full |> filter(country == countries[3]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = bio12 * `nitrogen_0-5cm_mean`) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Annual precipitation (mm) * Nitrogen (cg/kg)", y = "DOY", title = countries[3]) ``` ### Altitude - Precipitation ```{r} #| fig-width: 6 #| fig-height: 4.5 #| column: screen #| layout-nrow: 1 df_full |> filter(country == countries[1]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = alt * bio12) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Altitude * Annual precipitation (mm) ", y = "DOY", title = countries[1]) df_full |> filter(country == countries[2]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = alt * bio12) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Altitude * Annual precipitation (mm) ", y = "DOY", title = countries[2]) df_full |> filter(country == countries[3]) |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> mutate(var_inter = alt * bio12) |> ggplot(aes(x = var_inter, y = bloom_doy)) + geom_point() + scale_x_log10() + geom_smooth( method = "gam", formula = y ~ ns(x, df = 2), color = "red", size = 0.5 ) + labs(x = "Altitude * Annual precipitation (mm) ", y = "DOY", title = countries[3]) ``` ::: # Principal Component Analysis (PCA) - En estos gráficos aplico la función distinct() por latitud y longitud porque la misma coordenada (así sean en años diferentes) tiene el mismo dato de bioclimáticas, suelo y coberturas, ya que esta información no posee temporalidad. Por esa razón resumo con la **mediana del `bloom_doy`** (ya que esta variable sí cambia con el tiempo). - Aprovechando que existe alta correlación entre las variables predictoras, ejecuto análisis de componentes principales para reducir la dimensionalidad e intentar ver si en un espacio más pequeño de variables se exhibe algún patrón de comportamiento relacionado con la variable respuesta. - Para facilitar las visualizaciones discretizo la variable respuesta en los siguientes niveles (ordinales): - **`Q1`:** valores de `bloom_doy` inferiores al valor del cuartil 1: $DOY < Q1$. - **`Q2`:** valores de `bloom_doy` mayores o iguales al valor del cuartil 1 y menores al valor del cuartil 2: $Q1 \leq DOY < Q2$. - **`Q3`:** valores de `bloom_doy` mayores o iguales al valor del cuartil 2 y menores al valor del cuartil 3: $Q2 \leq DOY < Q3$. - **`Q4`:** valores de `bloom_doy` superiores o iguales al valor del cuartil 3: $DOY \geq Q3$. - La variable respuesta discretizada `doy_categ` la introduzco al PCA como variable suplementaría cualitativa. - Independientemente del país las floraciones que no superan el Q1 (`DOY < Q1`) y las que se tardan más del Q3 (`DOY >= Q3`) exhiben diferencias al graficarlas en los tres primeros componentes principales, quedando ubicadas en lugares diferentes. Es importante mencionar que el patrón es menos visible en Corea del Sur. Floraciones que están entre el el Q1 y Q3 tienden a ser similares. Esta diferenciación que se aprecia entre floraciones precoces y floraciones tardías podría constituirse en un indicativo de perfiles diferenciales en el nicho ecológico o entorno ambiental al cual están expuestos los cerezos. - En el caso de Japón se observa que el componente 1 está asociado de forma positiva con variables que recogen información de temperatura (bio1, bio5, bio10, bio11) y precipitación (bio12, bio13, bio16), indicando que valores mayores a cero en el componente 1 los podemos asociar con valores altos de temperatura y precipitación, justamente bajo estas condiciones observamos las floraciones rápidas (puntos azules), de otro lado, este componente está asociado de forma negativa con la latitud, la longitud, bio4 (estacionalidad de la temperatura), bio7 (rango anual de temperatura) y en menor medida con el nitrógeno, indicando que valores mayores a cero de este componente estarán relacionados no sólo con valores altos de temperatura y precipitación sino también con menor variación en la temperatura (bio4) y menor rango de temperatura anual (bio7), por tanto, el perfil descrito es el que aplica a las floraciones precoces, justo lo contrario ocurre con floraciones tardías, donde se esperan temperaturas y precipitaciones bajas, con mayores variaciones en la temperatura a través del año, también podríamos afirmar que plantas que están ubicadas más al sur de Japón tienden a presentar floraciones retardadas. Por la ubicación de los puntos también podríamos afirmar que las floraciones precoces en Japón están moduladas por cambios de temperatura, sin embargo, otros factores como la capacidad de intercambio catiónico (*cec*) y la densidad del carbono orgánico (*ocd*) podrían ser influyentes, ya que se observa que estas variables tiene correlación positiva con el componente 2, indicando que valores mayores a cero en el componente 2 están asociados con suelos altos en capacidad de intercambio catiónico y suelos con mayor densidad del carbono orgánico. - En Suiza también se evidencia la diferencia entre floraciones precoces y floraciones tardías, no obstante, la correlación de los componentes con las variables permite inferir patrones de comportamiento ligeramente diferentes a Japón y Corea del Sur. El componente uno tiene correlación positiva con variables que recogen información de temperatura (bio4, bio5, bio10, bio11) y exhibe correlación negativa con la altitud, el contenido de carbono orgánico del suelo (*soc*), el inventario de carbono orgánico del suelo (*ocs*) y en menor medida con variables de precipitación (bio12, bio16, bio17), estas correlaciones permiten inferior que floraciones precoses se dan en condiciones de alta temperatura (también la variación anual) y en sitios donde hay menor cantidad de lluvia anual, menor contenido de carbono orgánico del suelo y ubicadas en altitudes bajas. Los cerezos que exhiben floraciones tardías tienden a estar ubicados en altitudes elevadas, sitios con mayor cantidad de precipitación anual, mayores contenidos de carbono orgánico del suelo. - En Corea del Sur se observa un comportamiento similar al de Japón y Suiza en cuanto a que la temperatura es el factor determinante en los tiempos de floración. Sin embargo, el componente 1 está asociado de forma negativa con variables que recogen información de temperatura (bio1, bio6, bio9, bio11) y de forma positiva con variables que recogen información de la variación anual de la temperatura (bio4, bio7). Las floraciones tardías se ubican a la derecha del 0 en el componente 1, es decir, que estas plantas se espera estén en lugares con mayores variaciones en la temperatura a través del año y también se espera que las temperaturas de estos sitios sean más bajas. Basado en la asociación positiva que tiene el componente 3 con las precipitaciones (bio12, bio15, bio16) y observando la ubicación de los puntos verdes en el gráfico de CP1 vs CP3, es posible intuir que las floraciones demoradas también están asociadas con mayores precipitaciones. - Es importante mencionar que todas estas interpretaciones las hago desde la perspectiva exploratoria, siguiendo las premisas de [J. W. Tukey](https://web.stanford.edu/class/bios221/book/00-chap.html) como se muestra en el siguiente diagrama. La idea principal ha sido reconocer patrones que me permitan realizar un diagnóstico objetivo de cara a la fase de modelación. <center> ![](https://web.stanford.edu/class/bios221/book/imgs/iterativeparadigm.png){width=250} </center> ::: {.panel-tabset} ## Japan - Con cuatro componentes principales se retiene el 69.4% de la variabilidad total. - Con diez componentes principales se retiene el 89.2% de la variabilidad total. - Disminuir de 41 dimensiones a 4 perdiente aproximadamente el 30% de la variabilidad no está nada mal para propósitos de graficación, sin embargo, para propósitos de modelación, valdría la pena reducir de 41 variables a solo 10 componentes. ```{r} #| fig-width: 6.5 #| fig-height: 4.8 #| column: screen #| layout-nrow: 1 df_pca_japan1 <- df_full |> filter(country == countries[1]) value_q1 <- quantile(df_pca_japan1$bloom_doy, probs = 0.25, na.rm = TRUE) value_q2 <- quantile(df_pca_japan1$bloom_doy, probs = 0.50, na.rm = TRUE) value_q3 <- quantile(df_pca_japan1$bloom_doy, probs = 0.75, na.rm = TRUE) order_categ <- c("DOY < Q1", "Q1 <= DOY < Q2", "Q2 <= DOY < Q3", "DOY >= Q3") df_pca_japan2 <- df_pca_japan1 |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> select(-shrubs) |> mutate( doy_categ = case_when( median_doy < value_q1 ~ "DOY < Q1", median_doy >= value_q1 & median_doy < value_q2 ~ "Q1 <= DOY < Q2", median_doy >= value_q2 & median_doy < value_q3 ~ "Q2 <= DOY < Q3", median_doy >= value_q3 ~ "DOY >= Q3", ), doy_categ = factor(doy_categ, levels = order_categ) ) |> select(-c(location, country, bloom_date, bloom_doy, median_doy)) |> relocate(doy_categ, everything()) pca_japan <- PCA(X = df_pca_japan2, graph = FALSE, quali.sup = 1) df_pca_japan2$pc1 <- pca_japan$ind$coord[, 1] df_pca_japan2$pc2 <- pca_japan$ind$coord[, 2] df_pca_japan2$pc3 <- pca_japan$ind$coord[, 3] eigen_pc1 <- pca_japan$eig[, 2][1] |> round(digits = 1) eigen_pc2 <- pca_japan$eig[, 2][2] |> round(digits = 1) eigen_pc3 <- pca_japan$eig[, 2][3] |> round(digits = 1) df_pca_japan2 |> ggplot(aes(x = pc1, y = pc2, color = doy_categ)) + geom_point(size = 2.5, alpha = 0.75, shape = 18) + geom_hline(yintercept = 0, lty = 2, color = "firebrick3") + geom_vline(xintercept = 0, lty = 2, color = "firebrick3") + labs(x = glue("Component 1 ({eigen_pc1}%)"), y = glue("Component 2 ({eigen_pc2}%)")) + scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) + labs(color = "", title = "Component 1 vs Component 2") df_pca_japan2 |> select(-c(doy_categ)) |> correlate(method = "pearson") |> filter(term %in% c("pc1", "pc2", "pc3")) |> select(-c(pc1, pc2, pc3)) |> pivot_longer(-term) |> mutate( name = str_replace_all( name, "wildareas-v3-1993-human-footprint", "human-footprint93" ), name = str_replace_all( name, "wildareas-v3-2009-human-footprint", "human-footprint09" ) ) |> ggplot(aes(x = reorder(name, value), y = term, fill = value)) + geom_tile() + scale_fill_gradient2( low = muted("red"), mid = "white", high = muted("dodgerblue2"), midpoint = 0 ) + labs(x = "", y = "Component", fill = "Correlation") + theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7.5), legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) df_pca_japan2 |> ggplot(aes(x = pc1, y = pc3, color = doy_categ)) + geom_point(size = 2.5, alpha = 0.75, shape = 18) + geom_hline(yintercept = 0, lty = 2, color = "firebrick3") + geom_vline(xintercept = 0, lty = 2, color = "firebrick3") + labs(x = glue("Component 1 ({eigen_pc1}%)"), y = glue("Component 3 ({eigen_pc3}%)")) + scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) + labs(color = "", title = "Component 1 vs Component 3") ``` ## Switzerland - Con cuatro componentes principales se retiene el 66.4% de la variabilidad total. - Con diez componentes principales se retiene el 85.5% de la variabilidad total. ```{r} #| fig-width: 6.5 #| fig-height: 4.8 #| column: screen #| layout-nrow: 1 df_pca_switzerland1 <- df_full |> filter(country == countries[2]) value_q1 <- quantile(df_pca_switzerland1$bloom_doy, probs = 0.25, na.rm = TRUE) value_q2 <- quantile(df_pca_switzerland1$bloom_doy, probs = 0.50, na.rm = TRUE) value_q3 <- quantile(df_pca_switzerland1$bloom_doy, probs = 0.75, na.rm = TRUE) order_categ <- c("DOY < Q1", "Q1 <= DOY < Q2", "Q2 <= DOY < Q3", "DOY >= Q3") df_pca_switzerland2 <- df_pca_switzerland1 |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> select(-shrubs) |> mutate( doy_categ = case_when( median_doy < value_q1 ~ "DOY < Q1", median_doy >= value_q1 & median_doy < value_q2 ~ "Q1 <= DOY < Q2", median_doy >= value_q2 & median_doy < value_q3 ~ "Q2 <= DOY < Q3", median_doy >= value_q3 ~ "DOY >= Q3", ), doy_categ = factor(doy_categ, levels = order_categ) ) |> select(-c(location, country, bloom_date, bloom_doy, median_doy)) |> relocate(doy_categ, everything()) pca_switzerland <- PCA(X = df_pca_switzerland2, graph = FALSE, quali.sup = 1) df_pca_switzerland2$pc1 <- pca_switzerland$ind$coord[, 1] df_pca_switzerland2$pc2 <- pca_switzerland$ind$coord[, 2] df_pca_switzerland2$pc3 <- pca_switzerland$ind$coord[, 3] eigen_pc1 <- pca_switzerland$eig[, 2][1] |> round(digits = 1) eigen_pc2 <- pca_switzerland$eig[, 2][2] |> round(digits = 1) eigen_pc3 <- pca_switzerland$eig[, 2][3] |> round(digits = 1) df_pca_switzerland2 |> ggplot(aes(x = pc1, y = pc2, color = doy_categ)) + geom_point(size = 2.5, alpha = 0.75, shape = 18) + geom_hline(yintercept = 0, lty = 2, color = "firebrick3") + geom_vline(xintercept = 0, lty = 2, color = "firebrick3") + labs(x = glue("Component 1 ({eigen_pc1}%)"), y = glue("Component 2 ({eigen_pc2}%)")) + scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) + labs(color = "", title = "Component 1 vs Component 2") df_pca_switzerland2 |> select(-c(doy_categ)) |> correlate(method = "pearson") |> filter(term %in% c("pc1", "pc2", "pc3")) |> select(-c(pc1, pc2, pc3)) |> pivot_longer(-term) |> mutate( name = str_replace_all( name, "wildareas-v3-1993-human-footprint", "human-footprint93" ), name = str_replace_all( name, "wildareas-v3-2009-human-footprint", "human-footprint09" ) ) |> ggplot(aes(x = reorder(name, value), y = term, fill = value)) + geom_tile() + scale_fill_gradient2( low = muted("red"), mid = "white", high = muted("dodgerblue2"), midpoint = 0 ) + labs(x = "", y = "Component", fill = "Correlation") + theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7.5), legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) df_pca_switzerland2 |> ggplot(aes(x = pc1, y = pc3, color = doy_categ)) + geom_point(size = 2.5, alpha = 0.75, shape = 18) + geom_hline(yintercept = 0, lty = 2, color = "firebrick3") + geom_vline(xintercept = 0, lty = 2, color = "firebrick3") + labs(x = glue("Component 1 ({eigen_pc1}%)"), y = glue("Component 3 ({eigen_pc3}%)")) + scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) + labs(color = "", title = "Component 1 vs Component 3") ``` ## South Korea - Con cuatro componentes principales se retiene el 66.4% de la variabilidad total. - Con diez componentes principales se retiene el 89.7% de la variabilidad total. ```{r} #| fig-width: 6.5 #| fig-height: 4.8 #| column: screen #| layout-nrow: 1 df_pca_southk1 <- df_full |> filter(country == countries[3]) value_q1 <- quantile(df_pca_southk1$bloom_doy, probs = 0.25, na.rm = TRUE) value_q2 <- quantile(df_pca_southk1$bloom_doy, probs = 0.50, na.rm = TRUE) value_q3 <- quantile(df_pca_southk1$bloom_doy, probs = 0.75, na.rm = TRUE) order_categ <- c("DOY < Q1", "Q1 <= DOY < Q2", "Q2 <= DOY < Q3", "DOY >= Q3") df_pca_southk2 <- df_pca_southk1 |> group_by(lat, long) |> mutate(median_doy = median(bloom_doy, na.rm = TRUE)) |> ungroup() |> distinct(lat, long, .keep_all = TRUE) |> select(-shrubs) |> mutate( doy_categ = case_when( median_doy < value_q1 ~ "DOY < Q1", median_doy >= value_q1 & median_doy < value_q2 ~ "Q1 <= DOY < Q2", median_doy >= value_q2 & median_doy < value_q3 ~ "Q2 <= DOY < Q3", median_doy >= value_q3 ~ "DOY >= Q3", ), doy_categ = factor(doy_categ, levels = order_categ) ) |> select(-c(location, country, bloom_date, bloom_doy, median_doy)) |> relocate(doy_categ, everything()) pca_southk <- PCA(X = df_pca_southk2, graph = FALSE, quali.sup = 1) df_pca_southk2$pc1 <- pca_southk$ind$coord[, 1] df_pca_southk2$pc2 <- pca_southk$ind$coord[, 2] df_pca_southk2$pc3 <- pca_southk$ind$coord[, 3] eigen_pc1 <- pca_southk$eig[, 2][1] |> round(digits = 1) eigen_pc2 <- pca_southk$eig[, 2][2] |> round(digits = 1) eigen_pc3 <- pca_southk$eig[, 2][3] |> round(digits = 1) df_pca_southk2 |> ggplot(aes(x = pc1, y = pc2, color = doy_categ)) + geom_point(size = 2.5, alpha = 0.75, shape = 18) + geom_hline(yintercept = 0, lty = 2, color = "firebrick3") + geom_vline(xintercept = 0, lty = 2, color = "firebrick3") + labs(x = glue("Component 1 ({eigen_pc1}%)"), y = glue("Component 2 ({eigen_pc2}%)")) + scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) + labs(color = "", title = "Component 1 vs Component 2") df_pca_southk2 |> select(-c(doy_categ)) |> correlate(method = "pearson") |> filter(term %in% c("pc1", "pc2", "pc3")) |> select(-c(pc1, pc2, pc3)) |> pivot_longer(-term) |> mutate( name = str_replace_all( name, "wildareas-v3-1993-human-footprint", "human-footprint93" ), name = str_replace_all( name, "wildareas-v3-2009-human-footprint", "human-footprint09" ) ) |> ggplot(aes(x = reorder(name, value), y = term, fill = value)) + geom_tile() + scale_fill_gradient2( low = muted("red"), mid = "white", high = muted("dodgerblue2"), midpoint = 0 ) + labs(x = "", y = "Component", fill = "Correlation") + theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 7.5), legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm")) df_pca_southk2 |> ggplot(aes(x = pc1, y = pc3, color = doy_categ)) + geom_point(size = 2.5, alpha = 0.75, shape = 18) + geom_hline(yintercept = 0, lty = 2, color = "firebrick3") + geom_vline(xintercept = 0, lty = 2, color = "firebrick3") + labs(x = glue("Component 1 ({eigen_pc1}%)"), y = glue("Component 3 ({eigen_pc3}%)")) + scale_color_manual(values = c("dodgerblue3", "gray80", "gray80", "forestgreen")) + labs(color = "", title = "Component 1 vs Component 3") ``` :::